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We present numerical results for the local density of states in semiclassical Andreev 
billiards. We show that the energy gap near the Fermi energy develops in a chaotic billiard. 
Using the same method no gap is found in similar square and circular billiards. 
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I. INTRODUCTION 

The density of states in a normal metal in contact with a superconductor is affected by the superconductor, as a 
manifestation of the proximity effect or Andreev reflection. In early days of superconductivity, this effect has been 
observed in thin films of normal metal on a superconducting swbstrateJd It has been shown theoretically that for clean 
films .the spectrum of quasiparticle excitations remains gaplessa at Fermi energy, whereas a minigap develops for dirty 
films.13 The energy scale of this minigap is given by Ti/tN, ijv being a typical time spent by an electron in the normal 
metal before it gets to the superconductor.0 

Recent technological advances make it possible to study the effect in more complicated geometries, in diffusive 
metallic wiresflpjind in a two-dimensional electron gas where electron transport is almost ballisticO. Recent theoretical 
developmentsQ'tD suggested a new interpretation of old results. Bel The existence of a minigap has been related to chaotic 
character of the electron motion in the normal part of the system. It makes no qualitative difference whether the 
electron transport is diffusive, as in dirty films, or chaotically ballistic, as in clean billiardsa. The absence of the gap in 
the deGennes spectrum follows from the fact that a clean film is a specific case of a system with a separable geometry. 
In such a system the motion is not chaotic. 

This interpretation is rather difficult to comprehend in semiclassical terms. That was the motivation of the present 
research. The problem is as follows. The electron motion becomes truly chaotic (ergodic) only for trajectories that 
are very long in comparison with the system size. As we show in detail below, long electron trajectories correspond 
to Andreev levels with energies very close to Fermi energy. Therefore, in contrast to the interpretation in question, 
one may argue that the presence of a minigap shows that there are no long trajectories in the system. Consequently, 
it may not be chaotic. 
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FIG. 1. Panels a to c show the three clean billiards investigated. The specific points looked at are indicated in the chaotic 
billiard. The four normal/superconducting interfaces are labeled in the circular billiard. Panel d shows the form of a quasi ID 
diffusive system investigated as well. 



To resolve this sophistry, we have calculated the density of states for several Andreev billiards, chaotic and non- 
chaotic, depicted in Fig. [l]a to c. Such billiards combine Andreevoand specular reflection boundaries^ FIG. [l]a shows 
the form of the chaotic billiard investigated. At the circular parts of the boundary specular reflection occurs, while 
at the straight linear parts Andreev reflection takes place. The outward concave shape of the circular parts make the 
system highly chaotic. The special form chosen can be considered as representative for any chaotic Andreev billiard. 
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In order to show the marked difference with related integrable Andreev billiards, calculations are done also for circular 
and square boundaries. These systems are depicted schematically in FIG. [ljt» and c. The three systems differ only as 
far as the specularly reflecting boundaries are concerned. Their Andreev reflecting boundaries are identical. FIG. 
shows the form of a quasi one-dimensional diffusive system investigated as well. 

To find the density of Andreev states, we solve the equations for the quasiclassical Green's function along each 
classical trajectory. The solution depends explicitly on the length L of the trajectory considered, and gives a set of 
energy eigenvalues ~ Hvf/L. Then, for each given point, we calculate numerically all possible classical trajectories 
and sum up their contributions to the density of states. 

To put our results shortly, for the chaotic system we do observe the formation of a minigap near the Fermi level. 
Long, truly chaotic trajectories appear to take an exponentially small fraction of the phase volume and therefore do 
not contribute to the resulting density of states. The relevant equations are given in section ||. Results are discussed 
in section HI. Conclusions are formulated in section IV. 



II. THE EQUATION FOR LOCAL DENSITY OF STATES 

An expression for the local density of states n(e, r) in the clean, ballistic systems to be considered, will be derived by 
taking the imaginary part of the local Green's function G(e + iS, v, r). The quasiclassical Green's function G(iuj n , v, r) 
to be used will be a solution of the matrix equatior£3U 

v- VG + i[H,G] =0. (1) 

The velocity v is taken at the Fermi surface, the matrix H is given by 



H 



iu n A(r) 
-A*(r) -iwn 



(2) 



u> n are the Matsubara frequencies and A(r) is the superconducting gap function, to be taken constant in the super- 
conducting regions and zero in the central normal part of the system. Eq. ([!]) can be solved analytically for any 
trajectory. By accounting for all trajectories going through a given point the complete local solution can be found. 
We will show this by first writing Eq. (uh in the following from 



v F -^G + i[H,G] = 0, 

by which the r dependence is represented by the length parameter I along a trajectory. 
Suppressing the vp dependence, the general solution can be written in the form 



(3) 



G s (iu>n, t) = ciA s + c 2 B s e~ 
in which the matrices A$ and Bg are given by 
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A s = 
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-A* 
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c n + VK 2 + |A| 2 ) 



(4) 



(5) 



(6) 



while Ds = Bg. This solution is most easily obtained in two steps. First Eq. (||) is solved for the normal system, 
taking A(r) = in the matrix H . One finds 



Gjv(i^n, t) — c±A N + c 5 B N e "f + c 6 D N e "f 



with matrices 



A 



N 



1 

-1 



and B 



N 



1 





(7) 



(8) 
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while Djsr = B N . Since Eq. (||) is a homogeneous equation, the solution (]?]) is complete apart from an overall constant, 
to be determined by the requirement, that the matrix An times that constant is the solution of the original equation 
for the Green's function of p-bulk system, still having a delta function at the right hand side. This will merely lead 
to the proper normalizations of the expression for the density of states to be derived below. 
In the second step the full matrix H is diagonalized by the unitary matrix 



A — 7(^^+|Ap) V(l + "7(^f+|Ap) 



(9) 



The correspondingly transformed Eq. (||) has the same form as the equation for the normal system, the decay length 
being replaced by v F / V^n + |A| 2 )- Consequently, the full solution Eq. (ji|) of Eq. (|J) now can be obtained by the 
unitary transformation UGnU^ of the matrix given by Eq. (0), and substituting the proper decay length. 

The solution of Eq. (||) for a trajectory in inhomogeneous systems as depicted in FIG. [j] is obtained as follows. 
Consider an arbitrary trajectory, hitting some point at one of the superconducting/normal interfaces, at which point 
an Andreev reflection takes place, and follow the trajectory inside the normal region by accounting for all specular 
reflections at the boundaries, until it hits a superconducting/normal interface again. If the length parameter £ is taken 
at the initial hitting point and equal to L at the second hit, L standing for the total length of the trajectory, the 
form of the solutions in the different regions is clear from the Eqs. (||) and (0). For I < and £ > L the solution (||) is 
to be used, Since the dimensions of the superconducting regions are supposed to be large, behaving effectively as bulk 
superconductors, the coefficient of the matrix As can be taken to be equal to 1. Further, for £ < the coefficient C3 
has to be put 0, while for £ > L the B$ term blows up and the corresponding coefficient has to be put 0. The normal 
region is supposed to have mesoscopic dimensions, and the full solution G^i^n, £) has to be used. The equation (|^) 
being a first order differential equation, the only requirement is continuity at the interfaces. By that one finds for C4 
the following expression 



Un lAPtanh^ 
Ci - V^ 2 n + |A| 2 ul + |A| 2 + + I A| 2 ) tanh ' 



Now everything is ready for the local density of states in the normal region. First of all, since we are after studying 
the development of a gap just above the Fermi energy and of a width much smaller than |A|, it is sufficient to focus 
on the coefficient C4 in the limit |o; n | <C |A|, so on 

C4 = tanh . (11) 

v F 

Secondly, only the left upper matrix element of GN(iu n ,£) is required for the density of states, so only the C4 term 
in Eq. (Q) contributes. After the substitution iuj n — ► e + iS, one finds the contribution of a trajectory of length L 
through a chosen point, to the density of states at that point, to be proportional to 

i- rx... 1 (e + iS)L (e + iS)L ^-^ T/ eL , 

Inn 9ft tanh- — = hm9 ! tan- — = tt > S( (n+i)r), (12) 

5^0 iv F s^o v F ^ v F 2 

n 

in which the summation runs over integer n values. For a given point in the interior region of a 2D system all trajectories 
through that point have to be taken into account. This leads to the following expression for the dimensionless local 
density of states z/(e,r), defined by 

*fcr) S =fe£>= r^TA^-in+'M (13) 



flN Jo " v f 

in which tin is the constant density of states of a 2D normal system, and <f> is the slope angle of a trajectory with length 
L((f>). For presenting the results it is most convenient to shift to the dimensionless energy variable r\ = ei sys tem/wF, 
in which L sys tem stands for the length dimension of the system. By that, and denoting the relative density of states, 
now depending on the variable "q, by the same symbol v 1 we end up at the expression 

v( V , r) = /" S( V -^^ - (n + |)tt). (14) 

JO n ^system 

Note that in this final expression the length of a trajectory is present in a relative way, and has become dimensionless 
as well. 
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III. RESULTS AND DISCUSSION 



Because Eq. ( |l4| ) contains the relative quantity L((j)) / L systcm only, the physical size of the system does not enter. 
However, in actual calculations a choice has to be made. We have chosen L S ystem to be equal to 6. In all systems 
depicted in FIG. [I] the origin lies in the center of the billiard. The corners of the square billiard then lie at (±3, ±3). 
The length of the S/N interfaces is chosen to be equal to 2. The specific points 1 to 5 to be looked at have the 
coordinates 1 = (-2.2, 0.1), 2 = (-1.7, 1.4), 3 = (-1.2, 1.1), 4 = (-0.5, 0.1) and 5 = (-0.5, 0.2). 

In FIG. H the density of states is shown for point 1 in the three billiards under consideration, for a </>-grid with 
d(j) = 7r/200. For comparison the results according to finer grids are shown in FIG. |H[ In this figure the scale along 
the vertical axis has been adapted, by which the high peaks are not displayed fully, but the average structure comes 
out more clearly. The peak heights have been given explicitly. For the chaotic billiard, both figures show a gap at 
the Fermi energy, which corresponds with rj = 0. The missing states in FIG. || near r/ = for the square billiard are 
apparently due to the grid, because no gap is seen if ten and hundred times finer meshes are used. The finer mesh 
leads to a histogram, which is hardly distinguishable from a d<f> — 7r/2000 picture. But because of security all other 
histograms to be displayed are obtained using a dej) = 7r/20000 mesh. The density of states Square (??) is quite similar 
for the different points, and no particular features show up, so we further concentrate on the circular and chaotic 
billiards. 

The gap for the circular billiard is certainly typical for the point chosen. This becomes clear if one realizes, that 
states near 77 = are due to trajectories, which are long compared to the system size. For the circular billiard such 
long trajectories contribute only if the line through a chosen point and the center of the circle hits a normal, specularly 
reflecting boundary. This requirement is not fulfilled for point 1. At the top of FIG. || a similar plot is given, but for 
point 2, for which such trajectories certainly contribute. Now no gap is seen for the circular billiard. The gap for the 
chaotic billiard is manifestly present in the middle of this figure, for point 2, and at the bottom, for point 3. The gap 
in the histogram for the circular billiard depends critically on the precise location of the point. This is shown in FIG. 
f|, displaying results for the points 4 an d 5. While point 4 does not support long trajectories, point 5 does. For the 
latter point the gap is closed, although it lies very near to point 4- We show the histogram for the chaotic billiard at 
one of these points only, because not much difference is seen for this latter system. 

For the chaotic system the gap is present for all points. It appears to be an intrinsic property of this system. 
Long trajectories are rare, irrespective of the point considered. While for the square and circular billiards trajectories 
with over 2000 times the system size are easily found, for the chaotic billiard it is hard to find trajectories longer 
than 20 times the system size. This is illustrated by TABLE |j} The data given are obtained as follows. First, for 
fifty points the longest trajectory was calculated, using for each point a <f> grid of d<fi — 7r/2000. The longest of the 
100000 trajectories considered this way was found for point (-1.2, 0.9) at the angle given in the first line of the table 
in column 1. In addition to its length of 71.8 its relative length is given. In the last three columns the number of 
specular reflections and the labels of the S/N interfaces of both ends of the trajectory, called exits, are given. The 
meaning of these labels is shown for the circular billiard, in the middle of FIG. [l] After that the angle was specified 
finer and finer, using 10 angle values on both sides of the angle considered. Progressively the angle giving longest of 
the 21 lengths calculated that way was picked out for further subdivision. At the end, at one angle a length of 19.08 
times the system size was found, but then the borderline of our (double) precision was reached. The sensitivity to 
the initial condition is illustrated by the results on the 5th and 8th line from the bottom, because they hold for the 
same angle. This angle was generated in a slightly different way in the subsequent subdivision. We conclude that 
long trajectories are rare, although, theoretically they are supposed to exist. For example, comparing the trajectories 
illustrated by the third and fourth line, there must be an angle 4>o in between, at which the shift occurs from the first 
exit to the second one. This critical angle would support an infinitely long trajectory. 

To give some qualitative estimations, we consider the contribution of the trajectories nearing that critical angle 
(pa. The length at (f) — > <fio can be estimated as L{<f>) ~ — L systom log 2 (|</> — <jf>o | ) - This estimation follows from the fact 
that each bounce between concave boundaries approximately doubles the deviation angle. Using the relation ( |l4| ) we 
estimate v ~ 2~^. This shows that the density of states is exponentially suppressed at small r\. With our numerical 
method, having a finite grid 8<f>, we can only access energies 77 > — log 2 (<50). The states with smaller energies are not 
seen. The good convergence of our numerical data, even at relatively small 77, proves that the gap develops rather 
quickly, giving rise to an abrupt change of the density of states. 

It is interesting to note that quantummechanical effects can also be estimated in this way. Due to diffraction of 
the electron waves in a billiard geometry, the best angle resolution is limited by 8(f)n = 1/ ^k-p L systcm , ftp being 
the electron wave vector at Fermi energy. This implies that no Andreev states exists below the threshold energy 
77 ~ l/log 2 (fc F i sys tcm)- 

Although it is not the primary aim of the present paper, it is interesting to discuss the structure seen in the 
histograms. The peaks are certainly due to the fact, that special classes of trajectories contribute more than an 
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TABLE I. Illustration of the fact, that long trajectories are rare in the chaotic billiard. From top to bottom increasing 
lengths L(<f>) are found for an increasingly precisely chosen angle (j>. The data are for a point somewhat below point 3, with 
coordinates (-1.2, 0.9). 



in radians 



L(J>) 



L(<f>)/L sys 



No. reflections 



Exitl 



Exit2 



-0.92677083000000 



71.8 



11.96 



19 



-0.92677083600000 
-0.92677083700000 
-0.92677083800000 



75.5 
92.4 
70.8 



12.59 
15.40 
11.80 



21 

25 
20 



-0.92677083690000 
-0.92677083700000 
-0.92677083710000 



82.5 
92.4 
77.3 



13.76 
15.40 
12.88 



23 
25 
21 



-0.92677083700000 
-0.92677083701000 
-0.92677083702000 



92.4 
96.5 
77.8 



15.40 
16.09 
12.97 



25 
27 
22 



-0.92677083700000 
-0.92677083700100 
-0.92677083700200 



92.4 
110.3 
84.2 



15.40 
18.39 
14.04 



25 
30 
24 



-0.92677083700090 
-0.92677083700100 
-0.92677083700110 



88.9 
107.2 
91.0 



14.82 
17.86 
15.16 



25 
30 
25 



-0.92677083700096 
-0.92677083700097 
-0.92677083700098 



98.0 
114.5 
94.2 



16.33 
19.08 
15.69 



27 
32 
26 



average trajectory. Considering the L(ff) dependence of the local density of states expression (jl4|), for a given point 
and in an arbitrary direction this function will contain a term linear in <fi. But there are exceptions. The length L(4>) 
for trajectories through the points 1, 4 and 5 in FIG. |l| behaves quadratically in <f> around = 0, while for the latter 
two points in addition a quadratic behaviour in 4> — tt/2 around (f> — ir/2 holds. Then the argument of the delta 
function in Eq. ( fli] ) behaves quadratically in (f> around <ft = 0, giving rise to a square root singularity in the density 
of states. This singularity produces a series of equidistant peaks in the histogram with peak energies counted by the 
integer n. At 4> = the length L((f>) is equal to L S ystcm, so the lowest peak, for n = 0, is expected to be seen at 
rj = 7r/2 = 1.57. This peak is easily recognized in the FIGs. |^, |^ and ^. Another type of trajectory in a direction </>o, 
around which £(</>) behaves quadratically in cf> — cf>o, is, for example, the one through point 3 in FIG. [I] with </>o ss 7r/4, 
hitting the exits 1 and 4- Since the line from the origin to point 3 has a slope that is slightly smaller than -1, </)q is 
not equal to 7r/4. The corresponding L((j) ) = 5.23, leading to a peak at i] — 1.8, which is clearly present in the lower 
panel of FIG ra. The corresponding peak for point 2 lies too much to the right to be seen, namely at r\ = 2.3. We 
just point to another class of trajectories giving rise to quadratic behaviour, namely, for example, the one through 
point 3 and perpendicular to any outward concave circle. The corresponding line connects point 3 and the center of 
the circle. Although L(<fi) behaves quadratic as far as the contribution towards the circle is concerned only, possible 
linear contributions from the backward part of the trajectory will lead to a shift of the extreme (j>o. We calculated the 
corresponding lengths, and, just as an illustration, we mention the peaks for a few points. The trajectory through 
point 1 perpendicular to the circle centered at (3, 3) leads to a peak at 0.97, which is clearly seen in the lower panels 
of FIGs. U and [| The trajectories through point ,2 perpendicular to the circles centered at (3, 3) and (3, -3) lead to 
peaks at 0.97 and 0.49 in the middle panel of FIG. [|. 

Finally we want to point at even another type of extremal trajectory, namely a trajectory through a point and 
touching a circle. Consider the trajectory through point 2, in upward direction touching the circle centered at (-3, 3). 
The length has a minimum value at the touching angle 4>q, but the <j) dependence remains linear. Still an effect can be 
expected, because the coefficient of Acf> = cf> — cf>o for positive A</> can be different from the coefficient for negative Acf>. 
This leads to a step in the evaluation of the delta function in the expression for the local density of states, because 
6(bAcf>) = 6(Acf>)/\b\. In analyzing the different touching trajectories for the different points in the chaotic billiard, 
not all possible steps are clearly recognized, and others correspond to values of 77, which lie too high to be seen. We 
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just mention the trajectory through point 2 touching the circle around (3, 3), which is expected to give a step at 
r] = 0.77, and the trajectory through point 3 touching the circle around (-3, -3), corresponding to a step at 77 = 1.35. 
Both steps can be recognized in the middle and lower panel respectively of FIG |L 

We conclude by comparing the results for some point in the clean 2D chaotic billiard depicted by FIG. [l]a with the 
local density of states at a central point A and a point B closer to the S/N interface in the dirty ID Andreev billiard, 
depicted in FIG. [l]d. The results for the dirty system are obtained by numerical integration of the Usadel equations £3 
The curves A and B in FIG. ^ show the dimensionless density of states for the two points in the dirty ID system, 
while the histogram holds for point 2 of the clean system and reproduces the lower energy part of the middle panel of 
FIG. ||. Mind, that the scale of the gap energy E g is different in the ballistic and diffusive case: for a ballistic system 
Eg ~ Tivf I L system, whereas for a diffusive system it is strongly reduced, E g ~ hvp/ -wlLsysiem , I being the mean free 
path. In order to compare results, the curves and the histogram are rescaled to a common E g . 

Despite the gap occurs in both cases, it is seen that the behavior of the density of states is quite different for 
diffusive and ballistic cases. Consequently, this behavior is not universal and depends on details of the geometry and 
the scattering within the billiard. This fact strongly reduces the applicability of Random Matrix Theory methods, 
which are based on the assumption of universality of chaotic behavior. We note that the absence of universality can be 
understood from the fact that long and truly chaotic trajectories do not contribute to the density of states. Therefore 
it is determined by non-ergodic, non-universal trajectories. 



IV. CONCLUSIONS AND PROSPECTS 



Our results prove that a gap is formed in the density of states of a chaotic Andreev billiard even in the semiclassical 
limit. This is despite the fact that in the semiclassical limit Andreev states could have a very small energy being 
generated by very long trajectories. It turns out that the density of states of a chaotic billiard exhibits an abrupt drop 
at energies several times smaller than UF/Ajystem- Below this energy, the density of states is exponentially suppressed. 
We believe that quantum mechanical effects will lead to complete exhausting of the density of states at energies below 
VF/(isystcm log(fcFisystem))- Comparison of the density of states in the billiard and in the diffusive system clearly 
shows the absence of even approximate universality. 

For comparison a square and a circular Andreev billiard have been considered. In the square system long trajectories 
are always present, and no gap develops. Interestingly, although in roughly 70% of the volume of the circular system 
no gap is visible in the local density of states, in the remaining part of the volume there is also a gap developing. This 
possibly can be explained by the fact that the billiard is not ideally circular and may be slightly chaotic. 

All three billiards exhibit geometrically induced features in the density of states. Trajectories that traverse a billiard 
perpendicular to the superconducting boundaries, without reflection, generate a sequence of equidistant square root 
singularities. Trajectories that touch concave boundaries lead to steps in the density of states. 
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FIG. 2. From top to bottom the local density of states for the square, circular and chaotic billiards respectively, at point 1 
in FIG. El for the 0-grid indicated along the horizontal axis. 
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FIG. 3. The local density of states for the circular billiard at point 2, and for the chaotic billiard at the points 2 and 3 in 
FIG. 
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FIG. 4. The local density of states for the circular billiard at the points 4 and 5, and for the chaotic billiard at point 4. 




FIG. 5. The same as in FIG. |2| but for both a ten and hundred times finer (/>-grid for the square billiard, and for a hundred 
times finer grid for the circular and chaotic billiards. 
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FIG. 6. The local density of states at the points A and B in the ID dirty system depicted in FIG. |lp, compared with the 
histogram for point 2 in the 2D clean chaotic system depicted in FIG. ma. 
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